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ABSTRACT 

Observational and theoretical arguments suggest that the momentum carried in mass outflows from 
AGN can reach several times L/c, corresponding to outflow rates of hundreds of solar masses per year. 
Radiation pressure on resonant absorption lines alone may not be sufficient to provide this momentum 
deposition, and the transfer of reprocessed IR radiation in dusty nuclear gas has been postulated to 
provide the extra enhancement. The efficacy of this mechanism, however, will be sensitive to multi- 
dimensional effects such as the tendency for the reprocessed radiation to preferentially escape along 
sightlincs of lower column density. We use Monte Carlo radiative transfer calculations to determine 
the radiation force on dusty gas residing within approximately 30 parsecs from an accreting super- 
massive black hole. We calculate the net rate of momentum deposition in the surrounding gas and 
estimate the mass-loss rate in the resulting outflow as a function of solid angle for different black 
hole luminosities, sightline-averaged column densities, clumping parameters, and opening angles of 
the dusty gas. We find that these dust-driven winds carry momentum fluxes of 1-5 times L/c and 
correspond to mass-loss rates of 10-100 Mq per year for a 10® Mq black hole radiating at or near 
its Eddington limit. These results help to explain the origin of high velocity molecular and atomic 
outflows in local ULIRGs, and can inform numerical simulations of galaxy evolution including AGN 
feedback. 

Subject headings: black hole physics ~- galaxies: active ~ galaxies: kinematics and dynamics - galaxies: 
nuclei - radiative transfer - quasars: general 



1. INTRODUCTION 

1.1. Motivations from observations and theory 

The nature of the interaction between an accreting 
super-massive black hole (SMBH) and its host galaxy 
remains a challenging problem in the study of galaxy 
evolution. Numerical simulations reveal that gas can 
be drawn inward toward the nucleus by gravitational 
torques arising from a series of gravitational instabilities 



(Hopkins & Quataert 2010| ). This gas typically forms a 
dusty structure at small radii with a characteristic length 
scale of ^1-10 p arsecs which in some cases has been 



imaged directly (Jaffe et al. 2004 Raban et al 



•mm. 

led as 

a torus (iLawrence 1991 Antonucci 19931), but there is 



P henom enologicaily, this structure can b e mode 



an ongoing theoretical effort to provide a detailed, self- 
consistent explanation of its configuration and what sup- 
ports it. If a sufficiently strong poloidal magnetic field 
is present at the parsec scale, one possible explanation is 
that the dusty gas is launched as a hydrom agnetic wind 



( |Konigl fc Kartjel 1994| [Keating et al]|2012[ ). Heating of 
the ISM from stell ar feedback might s upport the dusty 
gas in a puffy disk f Hopkins et al.|20T2 ) . The disk might 
be simultaneously support ed by infrare d radiation pres- 
sure ( [Pier fc Krolik||1992||Krolik||2007| , or t he infrared 
radiation pressure ma y generate a tailed wind ( Dorodnit- 
syn et al.||2011 



2012) 
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Regardless of what supports the torus, gas continues 
to be drawn in to the black hole accretion disk at small 
radn (< 10"^'' cm), where it powers an active galactic 
nucleus (AGN). The radiation emitted from SMBH ac- 
cretion disks influences the dynamics of the torus itself, 
along with the dynamics of the host galaxy. This feed- 
back may act through a number of channels that include 
radiative heating (e.g. [Pi Mattco et al. 2005), jets (Silk| 



2005| [Croton et al.||20U5| |McNamara fc Nulsen||2007 



and winds dri ven by ramation pressure on resonan t ul- 
travio let lines ( Murray et al. 1995; Proga et al. 2000) and 



dust ((KoniKl fc Kartje 1994. ; Murray et al.^2005j Keating] 



et al. 2012p " Our challenge is to understand the com- 
bined ettect of all these modes of interaction. Improving 
our understanding of this connection will be crucial for 
answering questions about the growth of SMBHs, ob- 
servations of AGN, and the star formation histories in 
galaxies. 

Recent observations have begun to reveal the violent 
impact that AGN may have on their host galaxies. Ob- 
servations of obscured quasars such as SDSS J1356-I-1026 
have revealed outflows extending out to tens of kilo - 
parsecs from the galactic nucleus (Greene et al. 2012!). 
The estimated mechanical luminosity of these outflows 



(10 



,44-45 



ergs s ^) is too large to be explained by the 



inferred star formation activity. Other obscured quasars 
possess more massive outflows, with mass-loss rates of 
hundreds of solar masses per year ( |Moe et al. 



Dunn et al. 



2009 
loca 



2010|. Meanwhile, observations o: 
ultra-iummous infrared galaxies (ULIRGs) have led to 
the discovery of outflows with velocities th at are corre 



lated with the AGN bolometric luminosity ( Sturm et al. 
2011 ). These outflows also have mass-loss rates equal to 
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several times the star formation rate and in some cases 
exceeding 1000 solar masses per year, depleting the gas 
on timescales as short as 10® years. Adding to our picture 
are studies of post-starburst galaxies, exhibiting outflows 
with median velocity of approximately 1000 km s^^, sug- 
gesting t hat past AGN activity played a role in launching 
the gas dTremonti et al.||2007| . 



simulations of AGN feedback ( 


Ciotti et al.|2010 DeBuhr 


et al.|2011 


Debuhr et al.|201] 


J) that account tor deposi- 



radiation, including a combination of heating by X-rays 
and photoionizations, radiation pressure at the kilopar- 
sec scale, and winds driven from within a radius of less 
than 100 parsecs. Taken together, these effects ca n help 
to explain both the Mbh - o" relation (l^rr arcsc fc Mer- 
|ritt 2000;' Gebhardt et al.|2000|[^iinaiile'et al |2002^ ^ 
The existence ot galactic outflows observed at speeds of 
thousands of km s~^. The results, particularly those of 



Debuhr et al. (2011 1, also suggest that line-driven winds 
may be insutticient to drive observed outflows, and that 
a large amount of momentum (> 3L/c) may need to be 
deposited via absorption by dust grains during the pe- 
riod when the SMBH is optically thick to both ultraviolet 
and far-infrared radiation, t he time when most black hole 
growt h is believed to occur ( Fabian|[l999 Hopkins et al. 



20051. 

A large uncertainty in the numerical calculations ref- 
erenced above is the amount of radiative momentum de- 
posited within the central unresolved radius. The ve- 
locity and mass-loss rate of the resulting wind depend 
sensitively on this coupling. Moreover, in those studies 
the momentum was deposited in a spherically symmetric 
manner. In reality, multidimensional effects, such as the 
tendency for radiation to escape out the rarefied, polar 
regions of the gas distribution, will be crucial. These 
effe cts have been consid ered by several previous stud- 
ies. Pier & Krolik ( 1992 1 computed the radiation forces 
exerted on a toru s mode ll ed as a constant density cylin- 
drical shell, and |Krolik| (|2007|) extended that work to 
account for a more self-consistent rearrangement of the 
gas under the influence of the radiation. A radiation- 
hydrodynamics study that linked the effects of Gompton 
scattering and broad absorption line winds at the parsec 
scale with inflow processes on galactic (kilop arsec) scales 
in two spatial dimensions was undertaken by [Novak et al. 



(2011 1, and this was extended in ord er to capture the ra - 



diative transfer through dusty gas in Novak et al. (2012 ). 
Our study extends this work further by performing three 
dimensional Monte Garlo radiative transfer calculations 
for dusty gas, including both smooth and clumpy gas 
distributions, and by integrating the force on columns of 
gas in order to quantify the mass outflow rate from AGN 
radiating at high luminosity. 

The momentum flux in radiation from a SMBH accre- 
tion disk with luminosity L is L/c. Generally L will not 
exceed LEdd, the Eddington luminosity set by the elec- 
tron scattering (Thomson) opacity. Dust will contribute 
to the opacity seen by the radiation at large radii, but 
only at distances greater than the radius r^uh at which its 
temperature drops below the sublimation temperature 
Tsub ~ 1400 K. Although the subhmation temperature 
varies for each grain depending on its composition and 
its size, we choose to adopt the simpliflcation of assigning 



a uniform sublimation temperature to all the dust in our 
calculations. The sublimation radius may be estimated 
as 



L 



^sub 



47r(TsBT^ub 



0.62 pc 



L 



IQ46 ergs" 



1/2 



/ Tsub 

V 1400 K 



(1) 



When the gas distribution surrounding the SMBH is not 
isotropic, rgub may vary with angle. Within this radius, 
electron scattering dominates the opacity, and the usual 
Eddington limit applies. 

Once the intrinsic photons from the accretion disk en- 
counter dust in the surrounding gas, they are absorbed 
and the energy is re-emitted at infrared wavelengths. If 
the gas is also optically thick to the infrared, then the 
re-emitted radiation will continue to be absorbed and 
re-emitted in a random-walk pattern until it exits the 
optically thick region. Along the way, momentum will 
be imparted by the photons to the gas multiple times. 
In this scenario it is possible for the radiation to transfer 
momentum to the gas at a rate that exceeds L^dd/c For 
a spherically symmetric problem, this "boost" factor to 
the infrared radiation force is exactly the infrared opti- 
cal depth of the gas, which can be shown as follows: In 
steady-state, when radiative equilibrium holds and the 
luminosity as a function of radius is constant, we may 
compute the radiation force per volume /rad ^ 



/r 



L 



ad 



■ p{r) K{r) 



47r c 

The total outward force exerted by the radiation is 
L 



(2) 



/rad dV ^ 47r 



47r c 



OO 

p{r) K{r) dr ~ t— , (3) 
c 



where r is the radial optical depth for the infrared pho- 
tons. 

In a gas rich galactic nucleus with a column density 
of 10^^ cm~^, a mean mass per particle of 1.5 times the 
proton mass, and an infrared dust opacity of 10 cm^ per 
gram of gas, an initial guess for the optical depth would 
be approximately 250. There are two primary effects 
that will reduce the actual radiation force from such a 
high value. The first is the lack of spherical symmetry: 
a torus obscures only a fraction of the solid angle sur- 
rounding the accretion disk, and the presence of clumps 
and voids in the torus can increase the photon mean free 
path for certain sightlines. The second effect is dust sub- 
limation: dust will be absent from the innermost regions 
of the nucleus that contribute a substantial fraction to 
the gas column density, and the force integral can be 
well-approximated by setting its lower limit to rgub- 

To get a sense of the sort of momentum deposition 
rates that have been observed, consider the case of Mrk 
231. This system features an outflow of neutral gas with 
velocities in the range 360-900 km s~^ and a mass-loss 



rate estimated at 420 solar masses per year ( [Rupke fc 
Veilleux|2011 ) . The momentum flux in the outflow, esti- 
mated by multiplying the mass loss rate by the velocity, 
is between 2.6 to 6.5 times Ljc where L is measured to 
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be 1.1 X 10'*^ ergs s^^. The kinetic luminosity of the out- 
flow, on the other hand, is estimated at 7.3 x 10''^ ergs 
s~^, less than 1% of the bolometric AGN luminosity. 

Modeling the force from radiation pressure, and pre- 
dicting by what factor it exceeds L/c, becomes a difficult 
problem to tackle analytically in the absence of spherical 
symmetry, the presence of clumps, and with an account- 
ing for dust sublimation. For these reasons, we turn here 
to three-dimensional radiative transfer calculations using 
the wavelength-depen dent Monte Carl o radiative trans- 
fer code SEDONA (K asen et aL]|2006[ ). Given that the 
radiative diffusion time m these systems is shorter than 
the dynamical times, we restrict ourselves to steady-state 
configurations that do not include an explicit coupling to 
hydrodynamics. 

In section [2j we describe how we parametrize the gas 
configurations surrounding the black hole and how we 
treat the key physical processes in the radiative trans- 
fer. In section [3l we present our results for a series of 
calculations in which we vary the opening angle of the 
torus, the amount of gas present, and the accretion disk 
luminosity. We also examine how our dynamical conclu- 
sions are affected by accounting for a clumpy rather than 
smooth distribution of dust and gas. Finally, in section 
[4] we present our conclusions. 



2. METHODOLOGY 

2.1. Initial gas configuration - parameterized, smooth 

model 

Although the specific region we are studying is diffi- 
cult t o observe directly, gr avito-hydrodynamic simula- 
tions (Hopkins et al. 2012) provide information about 
its configuration before the effects of radiative feedback 
are felt. The gas and stars form a thick disk roughly in 
vertical hydrostatic equilibrium (our usage of the word 
"disk" throughout the remainder of this study refers to 
what is usually labelled as the torus and should not be 
confused with a reference to the black hole accretion disk, 
which is unresolved at our scales of interest ). The puffi- 
ness of the disk in the Hopkins et al. (2012) simulations 
is to some extent determined by the sub-grids turbulent 
velocity dispersion when strong stellar feedback in the 
ISM is included, but also by bending modes (firehose in- 
stabilities) driven by resolved velocities when less stellar 
feedback is included. While further accretion of the gas 
at this scale will rely on non-axisymmetric torques, we 
first adopt a simple axisym metric, hydrosta t ic disk model 
analogous to one used in 
parametrization captures t 



Hopkins et al.| (|2012| 

le key features o: 



This 
the gas con- 



figuration seen in the hydrodynamics simulations, but al 
lows us greater control over free parameters and removes 
unnecessary complications in our attempt to isolate the 
effects of the radiation. Such a parametrization also al- 
lows us to systematically introduce dumpiness into the 
gas for certain calculations (which, among other effe cts, 
breaks axisymmetry) , as will be described in section 2.2 
The vertical structure of the smooth disk model may 
be calculated by solving the equation of hydrostatic equi- 
librium in the normal (z) direction, assuming an isother- 
mal equation of state with effective sound speed Cs set by 
both the resolved and sub-grid velocity dispersion, along 
with any contribution from the thermal pressure of the 



(4) 



gas, 

Cs^ dp 

p dz dz ^ 
with solution 

p(i?,z)=p(i?,0)exp{c72[$(i?,0)-<i>(i?,z)]} . (5) 

Here $ denotes the gravitational potential, p denotes the 
density of the gas, and R is the cylindrical radius. If we 
assume that the gravitational potential is dominated by 
the mass of the central black hole Mbh at these scales, 
then the density distribution is 



p{R, z) = p(i?, 0) exp 




1 



- 1 



.(6) 

In the limit of small z/ R, this yields a Gaussian vertical 
structure. In this limit, the ratio of the squared sound 
speed to the squared Keplerian velocity Vc functions as 
the ratio of the disk scale height to the cylindrical radius, 
and for convenience we choose to define a parameter that 
makes this identification universal: 



hs 

1r 



GMbh 
R 



-1/2 



(7) 



Moderately large values of h/R (> 0.2) are sug- 
gested by the observed fraction of obscured versus un- 
obscured quasars, although gene rally this fraction cor - 
rclates strongly with luminosity ( Maiolino et al. |2007[). 
Meanwhile, typical values of hg/R found in Hopkins et al. 
(2012) range from 0.1 to 0.5. In this study we will con- 
sider hs/R in the range 0.1 to 0.35. 

Mid-IR interferometric observations find the mid-plane 
density may be well-fit with a power-law R^^ where 7 
lies within a range of approximately 0.4 to 1.4, with a 
tendency toward larger values for more luminous AGN 
(Kishimoto et al. 2011). This is also in agreement w ith 
the simulations presented in Hopkins et al. (2012) in 
which 7 w 1.5. We find that the results tor varying 7 
correlate strongly with the corresponding change in torus 
mass within the computational domain, and so we choose 
to capture variations in torus mass and column density 
by varying the normalization of the radial density profile 
(the parameter po described below), while fixing 7 to 1.5 
for all calculations presented in this paper. 

Converting to spherical polar coordinates r and 9, 
where 9 is taken to be zero along the z-axis, we obtain 



p(r, 9) = po 



r sin( 



exp[{hs/R)-^{sme-l)] 



(8) 

Here ro represents some inner cut-off radius where the 
density is po, and to prevent the radial column density 
from diverging we take p{r < ro) = pQ. 

One undesirable aspect of this model is that it leads to 
an accumulation of mass in the polar region of the disk, 
where sin0 is small. To correct for this, we allow the 
density profile to drop as a power law in the spherical 
radius r rather than in the cylindrical radius R. This 
amounts to dropping the factor of (sin^)^''', which is 
only significant far from the disk mid-plane. This leaves 



p{r, 9) = Po 



ro 



exp[{hs/RrHsm9-l)] , (9) 
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TABLE 1 

Fiducial parameters 











hs/R 


AfH ( cm-2 ) 


radial density i/^Edd 








power-law 7 




0.3 


3.4 X 1024 


1.5 1 


10« 



Note. — The first three parameters set the gas density distri- 
bution, while the last two set the relative strengths of the radiation 
pressure and gravity. The mean mass per partiele is always set to 
1.5 times the proton mass throughout this paper. Note that the 
column density presented in this table corresponds to integrating 
the gas density from large radii to a distance of 0.1 pc from the 
BH. The column density computed by integrating to the edge of 
the dust sublimation radius is 9.5 X 10 cm"'^ if the other fiducial 
parameters are fixed. 

w hich is quite similar to th e phe nomenoloRical models 
of iGranat o fc Danese ( 1994 ) and Efstathiou & Rowan- 
[Robinson ( [1995 ) that were used to explain the properties 
ol spectral energy d istributions obs e rved i n dusty AGN. 

The results from Hopkins et al. ( 2012[ ) indicate that 
hg/R does not change by more than a factor of order 
unity for all R. For simplicity, we take hg/R to be a 
constant for all R and allow it to vary as a free pa- 
rameter for different disk models. For all calculations 
in this study we assume a black hole mass Mbh of 10^ 
M0, and we parametrize the luminosity as a fraction 
of the electron-scattering Eddington luminosity for that 
mass. As mentioned above, we also vary po, which sets 
the sightline-averaged column density A^h- Unless stated 
otherwise, A^h corresponds to the column density inte- 
grated to a distance of tq = 0.1 parsecs from the central 
black hole. Also, unless A'h is being varied explicitly, 
po is set so that the sightline-averaged column density 
is 3.4 X 10^^ cm~^, with a mid-plane column density of 



1.0 X 10 cm~ ■ These values are consisten t with the 
calculations from Hopkins & Quataert (2010) of surface 
densities of 10^^ - 10^^ Mq kpc~^ for the central 10 par- 
secs surrounding the black hole. The fiducial parameters 
are summarized in Table [H 

For this smooth density model we use a two- 
dimensional grid with spherical polar (r, 9) coordinates, 
with logarithmic spacing in the radial coordinate and lin- 
ear spacing in the angular coordinate. Our resolution is 
192 radial zones and 64 zones for ranging from to 
7r/2, with an assumed symmetry for 9 — >■ tt — 9. The 
radial zones span radii ranging from tq — 0.1 pc to an 
outer radius rout = 10^" cm (w 32.4 pc). The 0.1 pc scale 
was chosen because it is a larger scale than the typical 
black hole accretion disk, but also small compared to the 
typical dust sublimation radius. We ignore all momen- 
tum deposition inside the 0.1 pc radius, and since nearly 
all the momentum deposition occurs at and beyond the 
sublimation radius, the exact choice of innermost radius 
has little effect on our results. Slices of the gas density 
for the mode l developed in this sect ion, along with a sim- 
ulation from Hopkins et al. (2012), are shown in Figure 

m — 

2.2. Initial gas configuration - clumpy models 

It has long been predicted on theoretical grounds that 
the dusty gas surrounding an accreting SMBH will not 
be smoothly distributed, but will instead form clumps 
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Fig. 1. — Top: An example of a slice through the smooth model 
density distribution with the fiducial parameters listed in Table [T] 
except A'^H = 1-0 X lO'^^ cm"^ (chosen to match the simulation in 
the bottom panel). Bottom: A density slice taken from a hydrody- 
namical simulation of gas accretion onto a central black hole (see 
Hopkins et al. (2 012| ). Note that the color map is truncated at 
n = 10'' cm~^, and the density in the model distribution continues 
to drop below this value. 



(Krolik & Begelman 1988). This prediction has been 



supported by observations such as the variabilit y of x- 



iti et al.||2002 ) as 


well c 


ma] Honig et al. 


2010 



ture exists concerning radiative transfer through clumpy 
torus models, with many prescriptions for generating 
clumpy density distributions from an underlying smooth 
density model and comparing the results to observations 
dNenkova et a l.||2002[ |Elitzur et al.||2004| |Honig et al 
|2006; Schart mami et ai |2008| [Stalevski et al.|2012[ [fley" 
mann fc Siebenmorgen||2012p . 

(Jur method for generating the cl umpy gas di s tribu- 



tion s most closely resembles th ose of Honig et al. ( 2006 1 



and Schartmann et al. (2008). We use a three dimen- 



sional grid and spherical-polar (r, 9, (p) coordinates, with 
logarithmic spacing in the radial coordinate and linear 
spacing in the angular coordinates. Our resolution is 128 



Radiation Pressure on the Dusty "Torus' 



5 





Fig. 2. — A clumpy gas distribution corresponding to the fiducial 
parameters in Table [T| and the clumping parameters for model 3 
in Table [2] Not pictured is the diffuse background gas. The white 
cube drawn at the center has a side length of 2 parsecs. 

radial zones, 96 9 zones for all 9 ranging between and 
77, and 192 cj) zones for all (j) ranging from to 2ti. The 
density of each clump in a given simulation is the same, 
and a preset number of clumps are placed on the grid. 
The clump positions are sampled from a probability dis- 
tribution derived fr om a smooth density distribution as 
described in section |2.1[ If two clumps overlap in posi- 
tion, their densities are added. Each clump's radius is set 
to a fixed number of grid zones in a given simulation, and 
the logarithmic radial spacing of the grid causes the size 
and optical depth of the clumps to grow with increas- 
ing distance from the SMBH. Overlaid on the clumps 
is a diffuse, smooth background gas distribution that is 
generat ed b y multiplying the density distribution from 
section 2.1 by 10~^. An example is pictured in Figure [2j 

For each clumpy gas distribution, the total mass in 
the computational domain was set equal to that of our 
fiducial smooth density model. The parameters varied 
between each clumpy gas distribution are the ratio of 
the clump diameter to the radial distance of the clump 
from the black hole {dc\/r), the gas density in the clump 
(rici), and the number of clumps in the simulation volume 
A^ci- The choices of Nc\ and dc\ set the average number of 
clumps per line of sight, which we compute by averaging 
over all [9, (j)) sightlines with a weighting to account for 
the solid angle subtended by each sightline. 

We have chosen four combinations of clumping param- 
eters to allow the average number of clumps per line of 
sight to t ake on values as low as 6 (in line with the re- 
sults from |Mor et al. ( 2009 )) to as large as 105 in order to 
demonstrate a transition to the smooth density models. 
These parameter choices are listed in Table [2] 

Ultimately we find that it is the volume filling fractions 
of the clumpy gas models that correlate most strongly 
with the integrated force exerted by the accretion radi- 



TABLE 2 

Clumping parameters 



Model # 


dcl/r 


Ticl ( cm 3 ) 




average number of 






clumps per l.o.s. 


1 


.24 


9.8 X 10" 


36864 


105 


2 


.12 


7.8 X 10^ 


36864 


25 


3 


.24 


7.8 X 10^ 


4608 


13 


4 


.49 


7.8 X 10^ 


576 


6.5 



Note. — See text for description of parameters 

ation. If we let /(r) denote the ratio of the volume oc- 
cupied by at least one clump to the total volume within 
a sphere of radius r centered on the black hole, then we 
find that it can be well approximated via broken power 
laws. For model 1, 



/(r) « 0.1 X 



0.05 X 



0.27 pc 



2.3 pc 



-0.33 



-1.5 



For models 2 through 4, 



/(r)«0.1x 



0.03 X 



0.19 pc 



-0.75 



1 pc 



-1.5 



for 0.1 pc < r < 2 pc , 

for r > 2 pc . 
(10) 

for 0.1 pc < r < 1 pc , 

for 7' > 1 pc . 
(11) 

Figure[3]shows the distribution of column density along 
randomly sampled sightlines for both our smooth and 
clumpy models. All column density values quoted in this 
study assume a mean mass per particle of 1.5 times the 
proton mass. 

Making the gas clumpy leads to a larger number of 
sightlines with lower column densities compared to the 
smooth gas distribution, and spreads out the peak on the 
higher end of the column density distribution. Both of 
these eff ects are more i n line wit h observa tional s urveys 



of AGN dRisaliti et al| l999 Ak ylas fc Ge organto poulos 
20091 [Malizia et al.||200'5( |LaMassa et al.||2009, ,Tr'eister 



et al. 2009). At the same time, our clumping prescrip 
tion tends to make the column density distribution bi 
modal, with a division between sightlines that intersect 
no clumps versus those that intersect at least one clump. 
This bi-modality, which is not present in the observa- 
tions, persists for all clumping parameters considered in 
this study, although it can be avoided if a larger fraction 
of the mass is allocated to the diffuse phase. 

2.3. Monte Carlo Radiative Transfer 

The Monte Carlo technique partitions the luminosity of 
the accreting black hole into equal-energy photon packets 
that probabilistically interact with the surrounding gas. 
The packets were transported in three dimensions for all 
calculations in this study. We improve our statistics by 
mapping the energy and momentum deposited by the 
packets into a two-dimensional array of zones - a photon 
that scatters at spherical coordinates {r,(j),9) is mapped 
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Fig. 3. — Top: Column density histogram (integrated for all gas 
from 0.1 pc to large radii) for a smooth density model with the 
fiducial parameters listed in Table [T] Middle: Clumpy gas column 
density distribution with same clumping parameters as those used 
in Figure [2| also integrated for all gas from 0.1 pc to large radii. 
Bottom: AU parameters are the same as the panel above, but this 
time the column density is integrated from the sublimation radius 
outward (i.e. these are the dusty gas columns). 



to position {r,0') where 6*' = 6* if < 6* < 7r/2 and 0' = 

TT- 61 if 7r/2 < 61 < TT. 

In Monte Carlo radiative transfer, the specific intensity 
of the radiation /(r, n, A) is constructed by counting the 
number of photon packets with wavelength A that enter 
into each grid zone at position r and with direction vector 
n in a given interval of time. Specifically, the radiation 
force per volume /rad at a given position is defined as 

/rad = - / pnxIndujdX. (12) 

To compute /rad in a given zone of our computational 
domain with volume AV over a time interval Ai, we per- 
form a sum a sum over all photon packets entering the 
zone. Each photon packet carries with it an energy Ep, 
a direction of travel rip, and a wavelength Ap. Associ- 
ated with that wavelength is an opacity K{Xp), measured 
per gram of gas, and which depends on whether dust is 
present at location r. The packet traverses a path of 
length Ar within a zone at position r. The force is then 

\ / p 

The radiative acceleration Orad is simply defined as 

/rad/P- 

Our calculations apply the stationarity approximation, 
in which we solve the steady-state radiative transfer 
problem for a fixed gas density distribution. This ap- 
proximation is justified if the radiative heating time scale 
and the radiative diffusion time scale are much shorter 
than the dynamical time scale. 

For a sound speed of 200 km s~^ and a characteris- 
tic length scale of 10 pc, the dynamical time is approx- 
imately 10^^ seconds. Meanwhile, the photon diffusion 
time through the disk never exceeds 10^^ seconds, and for 
many disk parameters the diffusion time is substantially 
shorter than that. The radiative heating time, estimated 
by dividing the thermal energy of the gas by the rate of 
radiative energy deposition, is 

^'''"^^ \ fj.mp J \pKcaT^\^J 

= ' (lofk) (loSlc) (lO cmVg) ' 

(14) 

which is also much shorter than the dynamical time. 

In this case, the condition of radiative equilibrium al- 
lows us to compute the dust temperatures by balancing 
radiative heating and cooling, 

47r J^pKi^hsiX) Bx{Tdust)dX 
^ J^ ^pKiihsWhdujdX . (15) 

For most calculations, the photons are emitted isotrop- 
ically at the edge of the 0.1 pc sphere surrounding the 
origin. The effect of anisotropic emission is treated in 
section [3^ We follow the photon propagation for time 
intervals of 5 x 10^ seconds, at which point we update 
the temperature of the dust in each grid zone. We treat 
dust as present everywhere where the dust temperature 
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is below 1400 Kelvin. The dust temperatures are up- 
dated until convergence is obtained at the one percent 
level, which typically takes fewer than 40 iterations if 
the initial dust temperature is set to 100 Kelvin in every 
zone. 

Finally, for estimating the dynamics of the gas based 
on the radiation pressure on the dust, we assume perfect 
hydrodynami cal coupling between the dust and the gas, 
as justified in Murray et al. (2005). 



og(rescaled a^^^) 



2.4. Intrinsic AGN spectrum 
We use the "intrinsic" (unr eddened) AGN spectr al en- 



ergy distribution described in Marconi et al. (2004). The 
majority of the spectral energy is found in the optical and 
near-UV and originates from the accretion disk, which re- 
sembles a 10^ Kelvin black body emitter. The spectrum 
also contains a sizable x-ray component. Intentionally 
absent from this spectrum is any infrared component, 
which we will calculate self-consistently based on the re- 
processing of the radiation by dust. 

2.5. Dust and electron interactions 

We use tabulated dust opacities and albedos based 
on Draine (20(]3a ) for wavel engths greater than 10 
Angstroms, and Draine (2003b) for shorter wavelengths, 
all corresponding to visual extinction ratio Ry = 3.1 
and assuming a fixed dust-to-gas mass ratio of 1/125. 
These values were interpolated between 48 reference 
wavelengths. In practice, the difference between scatter- 
ing and absorption is that for an absorption interaction, 
the wavelength of the re-emitted photon packet will be 
sampled from a probability distribution that depends on 
the dust's temperature, whereas the wavelength will re- 
main unchanged for a scattering interaction. For wave- 
lengths less than 100 Angstroms, we ignore scattering 
by dust since it will be almost entirely in the forward 
direction and hence will not lead to a net transfer of mo- 
mentum, although we still allow for absorption by dust. 

Electron scattering is only relevant for photons with 
wavelengths less than ^ 10 Angstroms, when the dust 
absorption cross section drops below that of the Thomp- 
son cross section, and when the photons are energetic 
enough to scatter equally well off of both bound and free 
electrons. We account for anisotropic, inelastic electron 
scattering in accordance with the Klein-Nishina formula. 

3. RESULTS 

3.1. Dust temperature and radiative acceleration 
dependence on smooth gas geometry 

Figure [4 shows slices of the equilibrium dust tempera- 
ture and the radiative acceleration vector field for disks 
of varying opening angles and with a smooth gas distri- 
bution. The color scheme is set so that all temperatures 
above the dust sublimation temperature appear as solid 
gray. Arrows representing the acceleration are plotted in 
zones where the dust is not sublimated and where the 
gas density exceeds 10~^^ g cm~^. 

We find that the dust sublimation region has an as- 
pherical, hour-glass shape. Sublimation extends to larger 
radii in the polar regions where the dusty gas is optically 
thin in the infrared. There, the dust absorbs ultraviolet 
radiation but emits in the infra red, forcing it to reach a 
higher temperature to maintain radiative equilibrium. 





I 




Fig. 4. — Arrows indicating the direction and strength of the 
radiative acceleration are plotted over slices of dust temperature. 
All parameters correspond to the fiducial values in Table n] except 
for opening angles which vary as indicated (while conserving mass 
in the calculation domain) . Regions in gray indicate where dust is 
sublimated (dust temperature that exceeds 1400 K). Acceleration 
arrows are present in zones where the dust is not sublimated and 
the gas density exceeds lO"'^^ g / cm"^. The arrow lengths are 
proportional to logio(10® X anet) where anet is in cgs units. 
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logtrescaled flux) 




R{pc) 

Fig. 5. — Arrows indicating the radiation fiux are plotted over 
gas density. All parameters in this calculation correspond to the 
fiducial values listed in Table ^ Arrows with significant deviation 
from the radial direction are colored black, while the boundary of 
the dust sublimation region is marked with black cells. The arrow 
lengths are proportional to logio(10~^^x net fiux (cgs)). Through 
a process of absorption of UV light and re-emission in the IR, at the 
inner wall of the dusty gas, fiux is channeled toward the poles in 
the outermost part of the dust sublimation region. The radiation 
travels radially in the dusty portion of the gas. 

Interestingly, nearly all the angular redistribution of 
the radiation occurs near the surface of the dust sub- 
limation region. Light from the central source initially 
travels isotropically to the inner edge of the dusty gas, 
and a large fraction of the photons are absorbed at the 
dust interface. When photons are re-emitted in the in- 
frared, many are sent back into the sublimation region. 
It is through this re-emission that the net flux becomes 
anisotropic at small radii. When infrared photons suc- 
ceed in penetrating deep into the dusty gas, they gen- 
erate a nearly radial radiative flux, as they would in a 
spherically symmetric problem (see Figure pi). 

Figure [6] displays how the radiative acceleration varies 
with radius and polar angle for the fiducial simulation. 
The behavior of the acceleration is quite different inside 
and outside the dust sublimation region - the presence 
of dust raises the opacity of the gas and therefore raises 
the radiative acceleration (as in equation 12). In a given 
solid angle the acceleration is highest just beyond the 
edge of the dust sublimation region, where ultraviolet 
and optical photons can push on optically thick, dusty 
gas. The acceleration rapidly drops as the radiation pen- 
etrates farther into the dusty gas and ultraviolet /optical 
light is converted into infrared, to which the dust is less 
opaque. For all solid angles, the acceleration settles to 
a constant ratio above gravity at sufficiently large ra- 
dius, indicating that the acceleration eventually becomes 
proportional to l/r^, further evidence that the infrared 
radiation diffuses primarily in the radial direction. In ad- 
dition to the radial dependence of the acceleration, there 
is an angular dependence that arises from the diversion 
of fiux from the mid-plane to the polar regions of the disk 
at the surface of the sublimation region. 

Slices of the net acceleration with gravitational accel- 
eration included are shown in Figure [7j In all cases the 
acceleration is primarily radial in direction, either out- 
ward or inward. Note that for opening angles hs/R < 0.3 



there is a critical polar angle below which radiation dom- 
inates over gravity and above which gravity dominates. 
In these cases infiow may persist in the equatorial region 
while gas is blown out at angles directed farther away 
from the mid-plane, potentially leading to a steady state 
outfiow. However, the radiative acceleration dominates 
over gravity everywhere when hs/R > 0.3 for this AGN 
luminosity and column density. 

For another perspective, in Figure M we plot the in- 
tegrated radiative acceleration for columns of g as as a 
function of polar angle (without gravitational accelera- 
tion included). We assume there are no forces in the 
tangential directions (i.e., each column is accelerated in- 
dependently), and that the radiation force is shared along 
the whole column as the inner gas pushes on outer gas. 
To compute this net acceleration, we first compute the 
integrated net force in each solid angle (including the 
effects of both gravity and radiation) , 



rad 



■ dr 



along with the mass in that solid angle. 



duj 



pr dr , 



(16) 



(17) 



where rgub denotes the edge of the dust sublimation re- 
gion for each value of 9. Then, the net integrated accel- 
eration is simply 



dFnct 




'dM^^' 


duj 


1 


duj 



(18) 



The value of Onct depends on t he choice of rout- How- 
ever, we will show in section 3.5 that the dependence of 



the rate of mass outfiow on rout is very small. 

From Figure [8] we see that as the opening angle of the 
parsec-scale disk becomes smaller, the radiative accelera- 
tion becomes more sharply divided between the optically 
thin and optically thick portions of the disk. This pri- 
marily reflects the sharper density gradients presen t for 
smaller opening angles. As we will show in section 3.2 
even though the radiative force is greater in the opti- 
cally thick portion of the disk, the force does not rise 
as quickly as the mass. This causes the acceleration to 
decrease toward the mid-plane. 

3.2. Enhancement of radiation force above L/c and the 
dependence on smooth gas geometry 

The ability of radiation to clear away ambient gas is 
enhanced by the fact that diffusing photons deposit their 
momentum multiple times as they random walk out- 
wards. We can quantify this effect in each solid angle 
by dividing the integrated force on the gas column in 
that solid angle by the radiative momentum per time per 
solid angle leaving the inner source. We call the resulting 
quantity Tcfr(0), which is computed as 



Teff(e) 







7 1 \ ^bh" 


duj 


1 


V 47r / c 



(19) 



and we extend the lower limit of the integral deflning 
dF-coAjditj from rgub to when computing this quantity. 
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Fig. 6. — Radiative acceleration in radius and solid angle. All parameters in this calculation correspond to the fiducial values listed 
in Table rn The acceleration is normalized at each radius by the gravitational acceleration GM^n/r^ , and the logarithm of that ratio is 
plotted. Thus, points with y-values above zero correspond to locations where the radiative acceleration exceeds gravity. The abrupt jump in 
acceleration occurs at the boundary of the dust sublimation region, where dust begins to contribute to the radiative opacity. As the radius 
increases beyond this boundary, the mean wavelength of the radiation transitions from the UV to the IR, rapidly lowering the radiative 
opacity in the process and reducing the radiative acceleration until eventually obeying an inverse square law dependence on radius. 



We may also compute an average value of this quantity 
averaged over all lines of sight, 



T eff 



1 

An 



Toff duj 



Tr/2 



Toff' sin 9 dO 



(20) 



Figure [9] summarizes our results for for various 
disk opening angles while holding the other parameters 
at their fiducial values as listed in Table [T] Increasing 
the opening angle boosts Toff (6*) for all 6, up to a max- 
imum value of approximately 5-6 for these parameters. 
In the polar region, this effect can be understood sim- 
ply in terms of the presence of more mass in that re- 
gion when the opening angle is larger. Meanwhile, even 
though there is less mass present in the equatorial region 
as the opening angle increases, the radiative flux in that 
region increases such that Tcs{Q) is able to increase there 
as well. 

If we calculate the radiation force on spherically dis- 
tributed gas with the same A^^h, we find that Tcff = 13. 
Thus, even though the effective radiation force exceeds 
L/c in Figure |9j the enhancement is not as large for a 
realistic disk geometry as it is in the spherically symmet- 
ric case. For the largest opening angle considered in this 
study {hs/R = 0.35), TcS is smaller than the spherically 
symmetric value by a factor of ~ 2. 

3.3. Results for Clumpy Gas 



the clumpy and smooth cases. This suggests that the 
smooth density distributions employed throughout most 
of this study provide accurate approximations to the be- 
havior of more realistic clumpy density distributions, al- 
though they should systematically overestimate the ra- 
diation force on the dusty gas by a factor of ~ 2 in the 
most extreme cases of clumping (fewest clumps per line 
of sight) considered here. 

Figure [TT] illustrates how the radiation force acts on 
portions of individual clumps. Note how the force re- 
mains radially directed even in the presence of clumps, 
and how clumps shadow gas behind them. 

3.4. Results for Anisotropic AGN Emission 

If the black hole accretion disk is aligned with the mid- 
plane of the gas present at the scale of our calculation, 
one might expect that there would be more flux emitted 
in the polar directions than i n the mid-pla ne direction. 
According to one prescription ( Netzer|1987 ), the emitted 
flux should obey 



^"01111110(1 oc cos 6* (1 -I- 2 cos 61) 



(21) 



where the first factor accounts for projected surface area 
and the second factor accounts of limb-darkening in an 
optically thick atmosphere. There is reason to doubt 
the validity of this model when relativistic effects are 
taken into account which tend to redirect radiation back 



Figure 10 shows how Toff (0) varies with the dumpi- 
ness of the gas. The shape of the momentum deposi- 
tion as a function of 9 appears generally the same for 



toward the mid-plane ( [Sun fc Malkan 1989 1. Moreover, it 
remains unclear whether the black hole accretion disk is 
aligned with the torus. We nevertheless choose to explore 
the scenario described by equation 21 in order to test 
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Fig. 7. — Arrows representing net acceleration (radiation + grav- 
ity) as a function of position. All parameters for this calculation 
correspond to the fiducial values listed in Table ^ except for open- 
ing angles which vary as indicated (while conserving mass in the 
calculation domain). Inward-directed arrows are colored black, and 
the arrow lengths are proportional to logio(10® X anot) where anct 
is in cgs units. For small opening angles, the gravitational acceler- 
ation exceeds that of the radiation in the equatorial region up to a 
critical angle above the mid-plane, and beyond this angle radiation 
dominates. As the opening angle increases, photons deposit more 
momentum in the dusty gas, and for sufficiently large polar angle 
the radiation force can exceed gravity in all directions. 
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Fig. 8. — Radiative acceleration as a function of polar angle (grav- 
ity not included). All parameters for this calculation correspond 
to the fiducial values listed in Table ^ except for opening angles 
which vary as indicated (while conserving mass in the calculation 
domain). The acceleration is lowest in the equatorial region, even 
though the force from radiation pressure is highest there, because 
the force does not rise as quickly as the mass as the polar angle 
increases. 
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Fig. 9. — Tcfr(9) for various opening angles (see equation |19| |. 
Parameters correspond to the fiducial values in Table ^ except 
for opening angles which vary as indicated (while conserving the 
total mass in the calculation domain). Teff ((?) is a measure of the 
radiative force on the gas column at a given polar angle, and it 
reaches its highest values in the equatorial region. 

the sensitivity of our resuhs to variations in the emission 
pattern of the accretion disk. 

Unlike the case of isotropic emission, we find that for 
anisotropic emission Teff(0) peaks at an intermediate an- 
gle < 7r/2. The peak arises because at small polar angles 
there is hardly any gas present to provide optical depth, 
whereas hardly any light penetrates into the dusty gas 
at large polar angles. For a calculation with our fiducial 
parameters, the ratio of Tog in the case of anisotropic 
emission versus Tes for the case of isotropic emission is 
0.72, indicating that photons tend to escape from the 
disk with fewer interactions when they are emitted in an 
anisotropic manner. This ratio will be even smaller for 
smaller disk opening angles. 
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Fig. 10. — t^b{6) for various clump densities n^i and clump sizes. 
Clump sizes are specified by the ratio of the clump diameter del to 
clump radial position r. In all simulations the number of clumps is 
varied such that the total mass of the gas in the simulations domain 
is held constant. The diffuse background makes up 1% of the mass 
in all simulations. The number of clumps per line of sight, listed 
in the legend, is computed by averaging over all (6, sightlines 
with a weighting to acount for the solid angle subtended by each 
sightline. 
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Fig. 11. — A slice of the magnitude of the radiation force for a 
clumpy gas distribution. A jump in the magnitude of the force is 
evident at the dust sublimation boundary. Arrows indicating the 
direction and magnitude of the radiation force are overlaid on one 
clump. The arrow lengths are proportional to logio{1.5 X 10^^) X 

/rad (cgs). 

If, instead of being diverted away form the torus plane 
as prescribed by equation |^ the radiation is beamed 
to^ward the plane, Teff would increase compared to the 
fiducial simulation. 

From this point on, we will only consider models with 
a smooth density distribution and isotropic central emis- 
sion. Nevertheless, it is important to keep in mind that 
the integrated force and M are likely to be modified due 
to the effects of gas clumping and anisotropic emission 
of radiation. 

3.5. Estimating the mass outflow rate 



We cannot determine precisely the dynamics of the 
gas without coupling the radiative transfer calculation to 
a hydrodynamic solver in a time-dependent calculation. 
However, we may apply an Eddington-type argument to 
approximate whether gas will be blown away in a given 
solid angle and to estimate the mass-loss rate. This ar- 
gument considers the gravitational and radiation forces 
but ignores centrifugal acceleration of the gas, viscous or 
gravitational torques, and shocks. The neglect of cen- 
trifugal support will not significantly alter the results in 
the cases when the radiation force on a column of gas 
is much less than or much greater than the correspond- 
ing force of gravity, but it will contribute to an under- 
estimation of the mass outflow rate in the intermediate 
range. 

Let t{6) denote the time taken to accelerate the gas 
in a given column with mass dMtot to a distance r^ut at 
constant acceleration anot(^)- Then 



m 
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a net 
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du 



dF„, 



(22) 



duj 



We define a differential mass outflow rate per solid an- 

yle dM /duj, 



dM 
duj 



(0) 



m \ 



dMg^ 
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(23) 



We can also deflne a mass outflow rate integrated over 
the entire volume (all of 9), 



M 



dM 
duj 



duj = 2 (27r) 



dM 
dw 



sin dO , (24) 



where we have taken advantage of the assumed symmetry 
for 6* — > TT — 0. Whenever dM/duj is less than zero for 
a particular value of 6, we do not add it to the total 
reported value for the total volume-integrated M, since 
we are only interested in the gas that get s blo wn away. 

Our gas density prescription (section 2.1) indicates 
that for any given polar angle, the density goes as r~^ . 
This allows us to compute dMgas/duj in terms of Tout and 
the sublimation radius rs^ihiO), 



dM,: 



duj 



cr /'(^-b) (^) 



r dr 
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In section 



3.6 



(25) 



we w ill prese nt values for dM /duj calcu- 



23 25 and the values of (iF„ct /d'-^ 



CTonte Carlo. For the rest of this 



lated using equations 
calculated using the 
section, we present a simple scaling argument to demon- 
strate that our estimates of the mass outflow rate depend 
only weakly on our choice of the outer radius rout (which 

is somewhat arbitrary). 

Our results from section [3T] indicate that we can think 
of the radiative acceleration as being divided into two 
parts: a spike in acceleration at the sublimation radius 
that arises from the absorption of ultraviolet and optical 
photons, and acceleration due to absorption of infrared 
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photons that goes as at large radii. Only the second 
type of acceleration is sensitive to our choice of Tout ■ We 
may approximate the infrared radiation force as 



dF„, 



du IR 

-7 r 



' ^ ' sub 



r dr 



' sub 
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(26) 



where arad(''sub) refers to the value of the radiative ac- 
celeration at the sublimation radius that provides the 
correct normalization for the inverse-square law acceler- 
ation at large ra dii. 

Using equation [23) dropping factors of order unity, and 
assuming Tout ^ ''subj we finally arrive at 
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(27) 



From the simulations performed in Hopkins et al.| 
(2012), 7 tends to fall between 1.5 and 2, and as already 
noted we have fixed 7 at 1.5 for all numerical calcula- 
tions in this study. We expect that the density profile 
will ultimately truncate at about 1 kpc. So, the ratio of 
the dM/du! due to absorption of infrared photons that 
we would calculate using rout of 1 kpc versus Tout of 32.4 
pc would be, for 7 = 1.5, only 2.4. For 7 = 2, dM/dut 
would be invariant with respect to choice of rout , and for 
7 = 2.5 the ratio would be 0.42. The fact that a signif- 
icant portion of the radiative acceleration in the Monte 
Carlo calculations comes from the spike near the dust 
sublimation radius further reduces the sensitivity of our 
results to our choice of rout ■ 

A final quantity that will be useful to us is the velocity 
of the gas in a solid angle f out (S) , 



t{9) = %/2 flnot ^'o 



(28) 



Once again focusing on the infrared acceleration at 
large radii and making the same approximations as we 
did for estimating dM /duj, we find 



{9) = j2,ri|-^ro 



^rad (^sub) 



GA/b 



^sub 
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(29) 



Therefore our calculations for the velocity of the gas will 
have a similarly weak dependence on rout as that of the 
mass-loss rate. 

3.6. Variation of mass outflow with opening angle 
Figure [12] shows the differential mass outfiow rate 
dM /dio calculated using equation 



23 



^for disks with vari- 
ous opening angles. The densities ot the innermost radial 
grid zones have been re-scaled to keep the total mass in 
the calculation domain constant in each case. 

The overall mass outflow rate declines with smaller 
hs/R due to the increased tunneling of radiation into 
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Fig. 12. — Differential mass outfiow rate dM /duj for various open- 
ing angles. Parameters correspond to the fiducial values in Table^ 
except for opening angles which vary as indicated (while conserving 
the total mass in the calculation domain). 

the low-density polar regions. For hs/R > 0.25 the dif- 
ferential mass outfiow rate peaks at 6 — 7r/2, while for 
hs/R < 0.25 the peak is at an intermediate polar angle. 
This is due to an interplay between the amount of mass 
available to be cleared away, its inertia, and the gravita- 
tional force acting upon it. Near the equator, however, 
the large amount of gas cannot be unbound by the radia- 
tive acceleration and so there is not outflow, even though 
the force due to radiation is strongest there. 

We emphasize that Figure [12] represents only a snap- 
shot in time of the mass outflow rate for an accret- 
ing SMBH. The evolution of dM/du: with time is not 
calculated here and requires a fully coupled radiation- 
hydrodynamics calculation. Depending on the resulting 
rearrangement of the gas, the long-term mass loss rate 
could conceivably be cither larger or smaller than the rate 
calculated here. One possible scenario is that an initially 
optically thick and puffy disk will blow away gas in the 
polar region. In the absence of a replenishing mechanism 
that operates on a timescale shorter than t, this might 
cause the disk to become thinner, reducing the tendency 
for radiation to blow out more gas. Alternatively, the re- 
moval of gas from above the near-midplane could change 
the geometry of the flux sufficiently so as to induce a 
stronger vertical component, drawing up more gas from 
the midplane and leading to yet more mass l oss. Yet an- 
other possibility, already suggested in section [XT] is that 
a steady-state infiow/outflow develops with a mass loss 
rate that hovers close to the instantaneous value com- 
puted here. 

3.7. Variation of mass outfiow with other parameters 



Figure 13 shows the differential mass outflow rate 
dM/du! for disks with varying column densities (the col- 
umn density is averaged over all lines of sight, and in- 
cludes both dusty and non-dusty gas). The variation in 
column density is directly proportional to variation in 
the total mass present in the calculation domain. 

As expected, more mass can be ejected when there is 
more mass present to begin with. However, we find the 

— 0.49 

scaling to be sub-linear: M oc for the range of col- 
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Fig. 13. — DifTerential mass outflow rate dM/dw for various 
sightline-averaged column densities (measured to 0.1 pc from the 
black hole). Parameters correspond to the fiducial values in Table 
[T] except for the average column densities which vary as indicated, 
xhe total mass in the calculation domain varies proportionally with 
the average column density. 



umn densities included in this study (power-law scaling 
relations for all the free parameters in the problem will 
be summarized in section 3.8). The slow growth of M 



with column density is due to the fact that at higher 
column densities, the radiative force on the gas in the 
densest portions of the disk does not rise as quickly as 
the mass present there, and so gravity becomes increas- 
ingly effective at lim iting the outflow rate. 

Finally, Figure 14 shows the differential mass outflow 
rate dM/duj for calculations with varying black hole lu- 
minosities. For higher AGN luminosities, not only is 
there a higher net force on a column at a given value 
of 6 for which the net force was already outward (posi- 
tive) , but also the net force becomes positive on columns 
at larger polar angles. For all other parameters held con- 
stant, there exists a critical luminosity at which the radi- 
ation force exceeds gravity for all polar angles, and all the 
gas would be blown away. For opening angle hg/R = 0.3 
and our fiducial mean column density 3.4 x 10^** cm^^, 
radial density power- law 7 = 1.5 and black hole mass 
Mbh — 10^ Mq, this critical luminosity is L/L^dd ~ 0.7. 
The existence of such a critical luminosity might help to 
explain the dearth of quasars observed to be radiating 
at the full value of their inferred Eddington limit, al- 
though the precise value of the limiting luminosity pre- 
sented here should be considered a rough estimate, and 
may vary with time as the gas rearranges itself following 
the initial outflow we have estimated. 



3.8. Summary settlings of integrated quantities 

The scalings of Tcff with the parameters of the problem, 
varied one at a time from the fiducial values listed, for 
a black hole with mass 10® Mq, can be summarized as 



to 
cu 
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Fig. 14. — Differential mass outflow rate dM/duj for various lu- 
minosities. Parameters correspond to the fiducial values in Table 
[1] except for the luminosities which vary as indicated. 

follows: 



TeS — 3.< 



hs/R 

0.3 



1.5 



0.49 



3.4 X 1024 cm-2 
L 

1.26 X 1046 ergs s^i 



-0.13 



(30) 



For these fits, six data points were used for hs/R span- 
ning 0.1 to 0.35, five data points were used for L spanning 
0.03 to 1. and four data points were used for A^h span- 
ning 102-^ to 3 X 1025 cm 2. All of these calculations 
used 7 = 1.5, and the results for changing 7 generally 
correlate with the results for the corresponding change 
in N-a- 

We may also present a summary scaling relation for 
the volume-integrated mass outflow rate M calculated 
over the same range of parameters: 



M = 144 Mq yr" 
TVh 



1 fhs/R 



0.3 



2.6 



3.4 X 1024 cm" 



0.62 



L 



1.26 X 1046 ergs s-i 



1.6 



(31) 



Note that A^h in these scaling relations corresponds to 
column densities integrated from 0.1 pc to large radii. If 
instead we use the column density integrated from the 
edge of the dust sublimation radius outward, then the 
fiducial column density becomes 9.5 x 102'^ cm~2^ the col- 
umn density power-law in equation [30] changes to 0.56, 
and the column density power-law in equation |3 1 1 ch anges 
to 0.71. Also note that the rates in equation |31| corre- 
spond to a radius of 32 parsecs from the ccntraTSMBH, 
and there is a weak dependence on radius (no stronger 
than r4/4 when 7 



3.5 



1.5; refer to section 
The fiducial value for the mass outflow rate of 144 Mq 
per year may seem surprisingly large. That value was 
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computed for a black hole radiating at its full Eddington 
luminosity, and at that luminosity the radiative accelera- 
tion beats out gravity at all solid angles. Therefore, there 
is reason to suspect that such a large outflow rate is not 
sustainable for many gas dynamical times at the parsec 
scale, as the amount of mass present and the opening an- 
gle of the disk readjust during the outflow. The quoted 
outflow rate also does not incorporate the effects of mak- 
ing the ga s dis tribution clumpy, and as was demonstrated 
in section |3.3[ this reduces the integrated force by a fac- 
tor of ^ 2 for significant clumping. Since the mass loss 
rate should roughly scale as the integrated force to the 
1/2 power (as argued in section 3.5 1, the mass outflow 



rate will be reduced approximately by a factor of \/2 in 
the case of significant clumping. On the other hand, the 
mass-loss rates calculated above correspond to the mass 
swept up within a radius of approximately 32 parsecs 
from the central black hole. Extra pol ating our results 
out to 1 kpc, as discussed in section |3l] might boost the 
outflow rates by roughly a factor of 2for 7 = 1.5. 

With those caveats in min d, th e i mpor tant conclusions 
to be drawn from equations [30] and [3l] are that the ra- 
diation force may reach several times (~ 3) L/c, and 
that the mass outflow rates can easily reach tens of so- 
lar masses per year for parameters close to our fiducial 
values. In the proper circumstances (A^h ^ ^O'^'^ cm~^, 
hs/R > 0.25, and L/L^dd ~ 1), the mass outflow rates 
can reach up to 100 A/0 per year. 

It is also interesting to note that the momentum 
enhancement and mass outflow rate depend relatively 
strongly on the AGN luminosity and disk opening angle 
hs/R. The fact that the dependence of the mass out- 
flow rate on luminosity is steeper than the dependence 
of the radiation force on luminosity may at first seem 
surprising . Th is scaling has its origins in an effect noted 
in section |3.7[ specifically that increasing the luminosity 
allows the radiation force to exceed gravity for a larger 
fraction of the solid angle, adding more mass to the out- 
flow than was present at lower luminosities. 

Finally, to drive home the point that momentum depo- 
sition, not heating, is responsible for the large computed 
outflow rates, we can estimate the corresponding kinetic 
luminosities. We use our estimate of the gas velocity as a 
function of solid angle Wout(^) (equation 28 1 to compute 
the fraction of the accretion luminosity L that goes into 
kinetic luminosity for the same parameter range used 
above: 



1 

L 

2(27r) 



1 dM 

2 da7 

''/^ 1 dM 

g 2 dlij 



{9)vl^,{9)doj 



{e)vl^M sine d9 



(32) 



Once again by varying each parameter one at a time 
with respect to the fiducial values, our results for a 10^ 



solar mass black hole can be summarized as 
1.9 



ek = 0.009 



K/R 
0.3 



3.4 X 1024 cm-2 



0.19 



L 



1.26 X 1046 ergs s"! 



1.8 



(33) 



If the column density is computed by integrating from 
the dust sublimation radius outward, the corresponding 
power-law in equation[33]barely changes at all, increasing 
to 0.21. 

The mass-weighted average velocity of the gas in the 
outflow will be approximately equal to TcffL/{cAI), al- 
though a more accurate value can be obtained by inte- 
grating Uout(^) weighted by dMgs^s/dO and only counting 
contributions from solid angles and radii for which gas 
can be blown out. For our fiducial parameters this av- 
erage velocity at 32 parsecs from of the computational 
domain is approximately 1000 km s^^. 

3.9. Comparison to Previous Results in the Literature 

It is important to compare and contrast the results 
of our calculation to previous studies that addressed 
the sa r ne phys i cal pr oblem ^ including |Pier fc Krolik| 



(19921, 



Krolikj ( |2007[ ) and [Dorodnitsyn et al.| ( |2"0TT 
I'he first of those studies included a calcula- 
tion of steady-state, multi-wavelength radiative transfer 
through a static dusty torus, but for a different torus 
geometry and slightly different boundary conditions for 
the radiative transfer than those used here. By modifying 
our de nsity distribution to match that of Pier & Krolik 



( 1992 1, and preventing dust from sublimating within the 
pre-determined boundaries of the torus, we find results 
that are in qualitative agreement with theirs, and quan- 
titatively the values for the components of the radiation 
force each agree to within 32% in the central region of 
the torus (our computed radiation forces are smaller). 
The primary difference between th e studies comes in the 
range of luminosities considered. |Pier fc Kr olik| ( |1992[ ) 
point out that for a large range of torus parameters, if 
L/L^dd ?J 0.1 the radiation force will overwhelm gravity 
and lead to an outflow, possibly halting accretion. This 
conclusion is consistent with our results, although we 
have sought to quantify the rate of mass outflow and have 
considered the possibility of simultaneous inflow and out- 
flow processes. 
The work of Krolik (2007) and Dorodnitsyn et al 



(2011 2012) generate self-consistent density distribu- 



tions for a torus that is dynamically influenced by radia- 
tion pressure, by assuming dynamical equilibrium in the 
first case and with a numerical radiation-hydrodynamics 
solver in the second and third. Again, these studies find 
that disruptive outfiows should occur when L/L^^^ > 
0.1 for Compton-thick torii, consistent with the present 
work. Further more, we find enco uraging quantitative 
agreement with Dorodnitsyn et al. (|2()12) when we use 
a gray radiative opacity of 10 cm^ per gram of gas to 
match theirs. For a 10^ Af© black hole with accretion lu- 
minosity of 10"*^ erg s~^ surrounded by a torus of mass 
5 X W^Mq enclosed within a radius of 3 parsecs, we find 
M using equations 23 and 24 to be 4.6 Mq per year, in 
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agreement with their value of 5 Mq per year. If we ex- 
tend the estimation of M to include gas out to a radius of 
32 parsecs, then our value of M rises to 27 Mq per year. 
Interestingly, when we perform the same calculation with 
wavelength-dependent radiative opacity, this value drops 
slightly to 25 Mq per year. This decrease is due in part 
to a larger region of dust sublimation in the wavelength 
dependent case. 

To summarize, our work is consistent with previous 
studies but seeks to quantify the mass lost in outflows 
driven by highly luminous (L/LEdd > 0.1) accretion 
events, accounting for mass driven away at large radii 
(> 3 parsecs). By integrating the momentum deposited 
in columns of gas, we can make a statement of how the 
outflows can affect the host galaxy. 

4. CONCLUSION 

We have calculated how radiation pressure from a lu- 
minous accretion disk around a SMBH drives a power- 
ful outflow of gas via continuum radiation pressure on 
dust at distances of 0.1-30 pc from the black hole. Using 
ambient gas conditions motivated by observational con- 
straints on nuclear obscuration in AGN [hg/R > 0.25, 
> 10^'' cm-2) we find that a 10* A/© SMBH radi- 
ating at Eddington can drive a wind with velocities of 
^ lOOO's of km s~^ and an instantaneous mass loss rate 
of 10-100 Mq per year (see equation [sT]) . For SMBHs 
with masses > IQ^Mq, the outflow rates could approach 

IOOOM0 per year. 

Radiative heating sublimates the dust out to distances 
of roughly 0.5 to 1 pc in the mid-plane, and radia- 
tion pressure drives away the gas and dust in the po- 
lar regions, leaving behind what may constitute the ob- 
served dusty torus. The wide-angle bipolarity of these 
outflows corresponds well t o observations of obscured 
quasars ([Greene et al. 2012) and Seyfert 2s (Crenshaw 
fc: Kraemer|2000 ). Although the radiative acceleration is 
greatest m the polar regions, the majority of the ejected 
mass comes from oblique angles where there is a more 
significant reservoir of gas. By contrast, gas in the equa- 
torial plane is more difficult to unbind because of its large 
inertia and large integrated gravitational attraction. 

The net momentum flux in the resulting outflow can 
exceed L/c by factors of up to 5 for the parameters 
studied, as infrared photons interact multiple times dur- 
ing their outward diffu sion . As recently demonstrated 
in the calculations of I Ciotti et al. | ([2010|; [Novak et al.| 
(2011 1; Debuhr et al. ( [2011 ), outflows with these proper- 
ties have a significant impact on gas in the surrounding 
host galaxy. Our results for the outflows match reason- 
ably well the observed outflows in lo cal ULIRGs such 
as Mrk 231 (?Rup ke fc Veilleux[ [20TT] ) . The mass-loss 
rates and kinetic luminosity tractions we calculate also 
provide a reasonable match to observations o f obscured 
quasars (Moe et al. 2009 Dunn et al. 20101, although 



our model does not provide a mechanism for launching 
large amounts of gas at the high velocities (> 20000 km 
s~^) observed in these systems at small radii. One pos- 
sibility is that these quasars are exhibiting both line and 
continuum radiation pressure driven outflows. 

We find that the net effect of the AGN radiation on 
the surrounding gas is a strong function of the opening 
angle of the accreting gas at the parsec-scale (the torus) . 



Increasing the opening angle allows more momentum to 
be deposited in all directions because the mass distribu- 
tion and emergent radiative flux become more isotropic. 
We also find a steep dependence of the outflow rate on 
the luminosity of the accretion disk, because at higher 
luminosities gas becomes unbound over a greater range 
of solid angles. This result is also in agreement with the 
observed anti-correlation between obscured AGN frac- 



tion and AGN luminosity 1 
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Keeping all of our parameters at the fiducial values 
listed in Table |l] but varying the luminosity, we find that 
outward radiative acceleration begins to exceed gravity 
at all angles once L/L^dd reaches a value of ~ 0.7. This 
value is subject to uncertainty given our approximate 
treatment of the gas dynamics, but it may nevertheless 
help to aid understanding of the relative dearth of broad- 
line quasars obs erved to be radiat ing at their full Edding- 
ton luminosity ( [Kelly et alT][2010[ ). 

The effects of dust sublimation play a crucial rule in de- 
termining the angular dependence of the radiative force 
on the torus. The redistribution of flux between polar 
angles takes place almost entirely in the region of gas 
in which dust has been sublimated, where infrared ra- 
diation is re-emitted from the edge of the dusty gas at 
angles deviating from the radial direction. Once photons 
penetrate into the dusty gas, they tend to diffuse radi- 
ally and deposit momentum almost entirely in the radial 
direction. 

All of the results presented above must be considered 
in light of the approximations and assumptions that we 
have used. In particular, while allowing for a range of 
torus scale heights and masses, we have focused on an ini- 
tial distribution of gas close to hydrostatic equilibrium, 
without accounting for the self-consistent dynamical re- 
arrangement of the torus as the outflow takes place, per- 
haps accompanied by inflow when possible. We have also 
discounted centrifugal support of the torus, which may 
lead to an underestimate of the mass loss rate when the 
integrated radiation force is close to that of gravity. 

Another effect we do not capture in these calculations 
is the potential for the outflowing gas to develop radiative 
Ray leigh- Taylor instabilities, which might provide more 
avenues of photon leakage and reduce the cou pling of the 
radiative mornenturn . to the gas as studied by [Krumholz 
& Thompson (2012). Concerning this last point, our 



clumpy gas simulations can provide a preliminary esti- 
mate of the extent to which the radiation force would be 
reduced in the presence of such instabilities. Also, due to 
the high opacity encountered by the direct ultraviolet ra- 
diation from the accretion disk, the direct radiation fleld 
plays a much more imp ortant role in launching th e gas in 
our calculation than in Krumholz fc Thompson ( 2012) . 
As demonstrated in Kuiper et al. "| 2012[ ), theradiative 
Rayleigh- Taylor instability can be suppressed when the 
direct radiation field is sufficiently strong. 

A fully coupled radiation-hydrodynamic calculation 
will be needed to fully understand the subsequent be- 
havior of the gas in time. Future work will focus on in- 
corporating the results of this study into hydrodynamic 
simulations of black hole accretion. In addition to the 
coupling to the hydrodynamics, more details pertinent 
to the radiative physics may be addressed in such calcu- 
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lations, including line absorption, anisotropic scattering 
off of dust, mctallicity gradients, and variations in the 
average dust-to-gas ratio. 
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